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Abstract 


A package  of  FORTRAN  subroutines  is  presented  for  the  solution 
of  nonsymmetric  or  symmetric  sparse  linear  systems  by  triangular 
decomposition.  Two  principal  aims  are  (1)  to  handle  matrices 
which  originally  fit  into  primary  core  storage  but  do  so  no  longer 
after  decomposition,  and  (2)  to  solve  a sequence  of  linear  systems 
all  of  which  have  the  same  sparsity  structure  by  generating --in 
secondary  storage- -a  record  of  the  decomposition  process  in  the 
form  of  an  integer  array.  Some  experimental  results  using  the 
package  are  included. 
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FURTHER  PROGRAMS  FOR  THE  SOLUTION  OF 
LARGE  SPARSE  SYSTEMS  OF  LINEAR  EQUATIONS1-* 

by 

21  21 
Charles  K.  Mesztenyi  ’ and  Werner  C.  Rheinboldt 

1.  Introduction 

In  a previous  report  [ 1 ] a package  of  FORTRAN  subroutines  was 
presented  for  the  solution  of  a linear  system 

(1)  Ax  = b 

based  on  triangular  decomposition  of  the  (symmetric  or  nonsymmetric) 
matrix  A.  The  underlying  data  structure  was  motivated  by  a more  general 
arc-graph  structure  discussed  in  [2  ]. 

The  programs  given  here  have  the  same  purpose  but  pursue  the  follow- 
ing two  different  aims: 

a.  To  handle  matrices  which  originally  fit  into  primary  core  storage 
but  do  so  no  longer  after  decomposition. 

b.  To  solve  a sequence  of  systems  (1)  all  of  which  have  the  same 
sparsity  structure.  This  case  arises,  for  example,  in  the  solu- 
tion of  nonlinear  systems  by  Newton's  method. 

The  first  aim  is  accomplished  during  decomposition  by  writing  the 
decomposed  part  of  the  matrix  into  secondary  storage  and  using  its  place 
for  newly  introduced  nonzero  elements.  In  order  to  meet  the  second  aim 

1"*This  work  was  supported  in  part  by  the  National  Science  Foundation  under 
Grant  GJ-35568JB  and  the  Office  of  Naval  Research  under  Grant  N00014-67-A 
-0239-0021  (NR- 044- 431) . 
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we  follow  an  idea  in  [ 3 ] and  generate--in  secondary  storage--a  record 
of  the  decomposition  process  in  the  form  of  an  integer  array.  This 
record  can  be  used  to  decompose  any  matrix  with  the  same  sparsity  struc- 
ture provided  there  are  no  round-off  problems. 


2.  Some  Background 


The  desired  triangular  decomposition  of  the  n * n,  nonsingular  matrix 
A has  the  form 


PAQ  = LU,  L = I + L , 


where  Lu  is  strictly  lower  triangular,  U upper  triangular,  and  the 
permutation  matrices  P,Q  define  the  pivoting  sequence.  The  decomposition 
is  accomplished  in  n steps,  such  that 


P^.  = (I+LpU.  + A.,  i = 0,1 n 


where  the  first  i rows  and  i columns  of  A.,  the  last  n-i  columns 

i’ 

of  L*?,  and  the  last  n-i  rows  of  IP  are  zero,  respectively.  Moreover, 

TO  0 t 

PP.L.  has  the  same  first  i columns  as  L and  U.Q.Q  the  same  first 

i rows  as  U.  This  latter  fact  allows  us  to  keep  L?  and  IL  in  second- 
ary storage. 

Let  t) (B)  denote  the  number  of  nonzero  elements  of  any  matrix  B. 

Then  the  storage  required  before  and  after  decomposition  is  of  the  order  of 
m0  = ^ m2  = h(U)  + respectively.  Furthermore,  m^  = max/n(Ap 

is  the  maximal  storage  needed  for  the  matrices  A^.  Clearly,  we  have 
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nig  2 m^  5 m?  and,  in  practice,  it  turns  out  that  n^-m^  veTT  much 
larger  than  m^-mp.  In  fact,  sometimes  we  found  m^  to  be  equal  to 
m()  (see  Section  5).  Hence  by  retaining  only  the  in  primary  storage 
we  require,  in  general,  only  little  more  storage  than  for  A itself. 

The  basic  storage  structure  allows  for  easy  modification  of  the  pivot 
ing  strategy.  In  fact,  in  the  nonsymmetric  case  the  pivot  selection  is 
handled  by  an  easily  replacable  subroutine.  We  use  here  the  well-known 
minimal  degree  algorithm.  If  is  the  set  of  nonzero  elements  of  A^, 
then  for  any  x ( S.  we  denote  by  R^(x)  and  CT(x)  the  subsets  of 
consisting  of  the  elements  in  the  same  row  and  column  as  x,  respectively. 
Now,  with  H^(x)  = R.  (x)  if  |R^(x)  | <)Ch(x)|  and  otherwise  E^(x)  = CN(x), 
the  set  of  potential  pivots  is  given  by 

(4)  S1?  = {x  * S. ; |a(x)|  > m-  max  |a(z)|}. 

1 1 z^Cx) 

Here  a(z)  is  the  value,  of  the  matrix  element  corresponding  to  z and 
M-  € [0,1],  a user-defined  parameter.  The  ith  pivot  is  then  the  element 
x * sj  for  which  ( |R^(x)  | -1)  ( |(T  (x)  | -1)  is  minimal.  Generally,  with  de- 
creasing p the  fill-in  decreases  while  the  round-off  influence  increases 

For  symmetric  A it  is  assumed  that  the  pivots  remain  on  the  main 

T 

diagonal  and  hence  that  Q = P . In  that  case  each  matrix  A^  is  again 
symmetric.  If  is  the  set  of  nonzero  diagonal  elements  of  A^,  then 
the  ith  pivot  is  the  element  x of  the  set 

(5)  D?  = (z  ( D.;  |a(z)  | 5 ii  max  |a(y)|) 

y®i 

for  which  the  number  of  nonzero  elements  in  its  row  is  minimal. 
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It  is  theoretically  possible  to  use  the  value  m-  = 0.  In  that  case, 

(4)  and  (5)  show  that  any  nonzero  element  of  or  D^,  respectively, 
is  a potential  pivot.  Then  the  pivot  selection  depends  only  on  the  sparsity 
structure  and  not  on  the  elements  of  the  matrix- -but,  of  course,  the  round- 
off influence  may  be  considerable. 


3.  Basic  Storage  Arrangements 

3.1  Matrices  in  Primary  Storage:  As  mentioned  before,  the  basic  storage 

structure  used  here  is  the  same  as  that  in  [ 1 ] . We  summarize  it  briefly. 

Set  N = {l,2,...,n}  and  let  S c N x N be  the  set  of  locations  corre- 
sponding to  the  nonzero  elements  of  a given  n x n matrix  A.  We  number 
the  elements  of  S from  n+1  to  n+m,  m = |S|,  that  is,  we  introduce  a 
bijective  mapping 

(6)  v:  S -*■  {n+1, . . . ,n+m)  . 


Now  define  two  integer  arrays  RY  and  CY  eacli  of  length  n+m  in 
which  the  relative  addresses  n+1,..., n+m  correspond  to  the  elements  of  S 
in  the  order  provided  by  v.  The  images  v(R^)  of  the  row  sets 


R^  = { (i, k)  i S;  some  k € N},  i ( N 


form  a partition  of  {n+1 n+m).  For  any  set  v(Rp  = {i^,...,i^}  we  link 

the  locations  i,^,^,...,^  into  a circular  list 


(7) 


Lx  = RY(i),  ij+1  - RY (ij ) , j = l,...,k-l,  i = R(ik) 


where  for  practical  reasons 
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Analogously,  we  proceed  with  the  images  v(C^)  of  the  column  sets 


Cj  - {(k,j)  f S;  some  k f N},  j ? N 

in  the  array  CY. 

In  order  to  store  the  associated  matrix  elements  a third  array  A is, 
of  course,  needed.  Thus,  for  example,  the  matrix 

0 1 0 3 

0 0 5 0 

-10  2 0 
0 -2  0 0 / 

may  be  stored  as  follows: 


loc 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

RY 

10 

8 

9 

7 

3 

1 

4 

2 

5 

6 

CY 

5 

7 

9 

10 

1 

2 

6 

3 

8 

4 

A 

* 

* 

A 

A 

-1 

1 

-2 

5 

2 

3 

We  shall  refer  to  RY  and  CY  as  the  sparsity  structure  arrays  and  to 
A as  the  coefficient  array. 

For  symmetric  A the  set  S only  needs  to  be  the  set  of  locations 
of  the  nonzero  elements  in  the  upper  (or  lower)  triangle  (including  the 
diagonal)  of  A.  Moreover,  we  assume  always  that  in  the  symmetric  case 
all  diagonal  elements  are  nonzero.  Then  the  first  n cells  of  the  sparsity 
structure  arrays  RY  and  CY  are  no  longer  needed  if  (6)  is  changed  to 
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v:  S -*■  {1,2,. ..,m),  v(i,i)  * i,i=l,...,n  . 


There  is  no  need  to  repeat  the  details;  the  resulting  data  arrangement 
should  be  self-evident  from  the  following  example: 


1 0-1  o\ 

0 2 0 -2 
-1  0 3 0, 

\ 0 -2  0 4 / 


loc 

1 

2 

3 

4 

5 

6 

RY 

5 

6 

3 

4 

1 

2 

CY 

1 

2 

5 

6 

3 

4 

A 

1 

2 

3 

4 

-1 

-2 

During  decomposition,  this  storage  arrangement  is  used  for  the  matrices 
A^,  i = 0,...,n.  When  a nonzero  element  of  remains  in  A^+j  its  posi- 
tion in  the  RY,CY,AY  arrays  is  maintained.  After  each  pivoting  step  the 
pivot  row  and  column  are  written  on  secondary  storage  and  the  corresponding 
cells  in  the  storage  arrays  are  freed,  that  is,  the  circular  linkages  con- 
taining these  elements  are  modified  appropriately.  The  resulting  free 
locations  are  reassigned  when  fill-in  occurs. 

3.2  Triangular  Matrices  in  Secondary  Storage:  The  programs  here  are 

written  for  use  with  a random-access  secondary  storage  device.  Some  infor- 
mation about  the  necessary  I/O  routines  is  provided  in  Section  4.1  below. 

In  the  nonsymmetric  case  the  triangular  matrices  L and  U are 
written  as  two  arrays  of  pairs  of  numbers.  The  L-array  contains  the  columns 
of  L in  consecutive  order  and  each  column  has  the  form 


* i - 1— 

c r 


i 


1 "'W 


where  ic  and  ir  are  the  column-  and  row- index,  respectively,  of  the 

pivot  (stored  negatively)  and  j-  represents  the  row  index  and  l.  • 

1 Ji’  c 

the  value  of  each  nonzero  element  in  the  particular  column.  Similarly  the 

U-array  contains  the  rows  of  U in  consecutive  order,  each  of  them  in  the 


J k,ui  i 
11  1r*^k 

'VPi 


Here  ic  is  the  column  index  of  the  pivot  and  p^  its  value,  while 
j.,u.  . denote  the  column  index  and  value,  respectively,  of  each  non- 

1 VJi 

zero  element  in  the  row.  The  entire  U-array  is  initialized  by  a dummy 
pair  -1,-1. 

For  the  backsubstitution  programs  it  is  assumed  that  the  L- array  is 
read  forward  and  the  U-array  backward. 

In  the  symmetric  case,  there  is,  of  course,  no  need  for  both  the  L- 
and  U-array.  Accordingly,  only  the  L-array  is  set  up  containing  the 
columns  of  L in  consecutive  order,  each  one  in  the  form 


-i -,d. 
d l. 


^k’^i  i 

-i,,d. 
d ij 
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Here  ij  is  the  index  of  the  pivot  (on  the  diagonal)  and  d^  its  value 
and  j.,l.  ■ have  the  previous  meaning.  The  header  at  the  beginning 

1 h,Xd 

and  end  of  each  column  is  needed,  since  during  backsubstitution  the  array 
is  read  once  forward  and  once  backward. 

3.3  Representation  of  the  Decomposition  Record:  As  mentioned  in  the 

introduction,  the  programs  can  generate  a record  of  the  decomposition  • 
process  for  later  use  with  any  other  matrix  of  the  same  sparsity  type. 
This  record  is  in  the  form  of  an  array  of  positive  integers  in  secondary 
storage.  For  each  pivoting  step  the  following  information  is  recorded: 
NonAytmeJyUc  Ccu>e.: 


kj. , Yy  m^, 


• • • 


l 


V 


i relative  location  of  the  pivot  in  the  RY,CY  arrays 

ic,  ir  the  column-  and  row- index  of  the  pivot,  respectively 
kc,  kr  the  number  of  nonzero  elements  in  the  pivot-column  and  pivot- 
row,  respectively 

xi*  -*i  re^attve  location  (in  CY)  and  row- index,  respectively,  of 
the  nonzero  elements  in  the  pivot  column 
y^,  m^  relative  location  (in  RY)  and  column- index , respectively,  of 
the  nonzero  elements  in  the  pivot  row 
relative  locations  of  the  elements  in  which  must  be 
modified  at  the  step,  t = k • k 

r»  c r 


1 


’T'  - 
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SymeX/Uc.  Ccu>&: 
i»  k>  y^, 


yk’  \ 


V \ 


l. 

J 


index  of  the  pivot  and  hence  also  its  relative  location  in  the 
RY , CY  arrays 

the  number  of  nonzero  elements  in  the  pivot  row.  Each  of  these 
elements  is  identified  by  a pair  (y^ ,nr ) as  in  the  nonsymmetric 
case 

the  relative  locations  of  the  elements  in  A.  to  be  modified, 

1 

t = k(k-l) 


4.  Description  of  the  Programs 

The  package  consists  of  four  groups  of  subroutines  with  names  INT,  BLD, 
DEC,  and  SLV;  in  addition,  there  is  a pivot  selection  routine  PVT01  for 
the  nonsymmetric  case  and  a set  of  I/O  routines  for  interface  with  the  random 
storage  device. 

The  INT  programs  initialize  the  storage  area  and  have  to  be  called  first. 
The  BLD  routines  establish  the  data  structure  described  in  Section  3.1  for 
the  given  matrix  A.  Then  the  DEC  routines  are  called  to  perform  the 
decomposition  of  the  matrix  and/or  to  generate  a record  of  the  decomposi- 
tion process.  Finally,  if  applicable,  the  SLV  routines  are  used  to  obtain 
the  solution  of  the  given  system  (1)  by  backsubstitution. 

In  general,  any  efficient  routine  for  building  up  the  basic  data 
structure  from  given  data  about  the  matrix  depends  strongly  on  the  details 
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of  the  files  used.  The  BLD  programs  presented  here  avoid  all  assumptions 
about  file  formats,  etc.,  by  establishing  the  data  scructure  one  matrix 
element  at  a time.  In  other  words,  the  chosen  BLD  routine  has  to  be  called 
once  for  each  nonzero  matrix  element.  For  many  practical  purposes  this 
may  be  inefficient.  The  routines  were  included  principally  for  the  sake 
of  completeness;  it  should  be  easy  to  rewrite  them  for  any  specific  appli- 
cation. 

The  names  of  all  subroutines  in  the  four  principal  groups  are  preceded 
by  the  letters  S or  N for  the  case  of  symmetric  or  nonsymmetric  matrices, 
respectively.  The  names  of  the  subroutines  in  the  INT,  BLD,  DEC  group 
are  ending  with  one  of  the  numerals  0,  1 or  01.  This  indicates  the  follow- 
ing alternatives: 

0 - Initialize  or  build  only  the  sparsity  structure  arrays  of  the  matrix 

or  generate  a record  of  the  decomposition  based  solely  on  the  sparsity 
structure.  These  routines  are  only  available  for  symmetric 
matrices;  for  nonsymmetric  matrices  it  is  not  an  advisable  approach 
since  the  resulting  round-off  error  could  be  severe. 

1 - Initialize  or  build  only  the  coefficient  array  for  the  matrix  elements, 

or  decompose  the  matrix  using  a previously  generated  record  of  a 
decomposition  for  matrices  with  the  same  sparsity  structure. 

01  - Initialize  or  build  both  the  sparsity  structure  arrays  and  the  coeffi- 
cient array,  or  decompose  the  given  matrix  and,  optionally,  generate 
a record  of  the  decomposition. 

The  pivot  selection  for  the  nonsymmetric  case  is  performed  by  the  rou- 
tine PVT01.  For  the  symmetric  case,  pivot  selection  is  incorporated  within 
the  routines  SDECO  and  SDEC01. 
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4.1  Catalog  of  Subroutines:  In  this  subsection  we  list  the  various  sub- 

routines of  the  package  together  with  their  calling  sequences  and  brief 
descriptions  of  their  purposes.  The  arguments  in  the  calling  sequences 
are  discussed  in  the  next  subsection. 

INT  - Routines 

SINTO(MD,RY,CY,ND) 

Initialize  the  sparsity  structure  arrays  of  a symmetric  matrix. 

SINT01  (MD,FD,RY,CY,A,AN,ND) 

Initialize  the  sparsity  structure  arrays  and  the  coefficient 
array  of  a symmetric  matrix. 

SINT1  (MD,FD,AN)  # 

Initialize  the  coefficient  array  of  a symmetric  matrix. 

NINT01  (MD,FD,RY,CY,AN,NDR,NDC) 

Initialize  the  sparsity  structure  arrays  and  the  coefficient 
array  of  a nonsymmetric  matrix. 

NIMT1  (MD,FD,AN) 

Initialize  the  coefficient  array  of  a nonsymmetric  matrix. 

BLD  - Routines 

All  routines  add  a matrix  element  with  value  V,  row  index  I,  and 
column  index  J to  the  structure.  Note  that  in  the  symmetric  case  only 
the  nonzero  elements  in  the  upper  (or  lower)  triangle  and  the  diagonal 
should  be  given. 
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SBLDO(I,J,MD,RY,CY,ND) 

Insert  element  (I,J)  into  the  sparsity  structure  arrays  of  a 
symmetric  matrix. 

SBLD01(I,J,V,MD,FD,RY,CY,A,AN,ND) 

Insert  element  (I,J)  into  the  sparsity  structure  arrays  of  a 
symmetric  matrix  and  a corresponding  value  V into  the  coeffi- 
cient array. 

SBLD1(I,J,V,MD,FD,A,AN) 

Associate  a value  V to  element  (I,J)  of  a symmetric  matrix. 

The  V-values  must  be  in  the  same  order  in  which  the  (I, J) -values 
were  read- in  during  prior  construction  of  the  corresponding 
sparsity  structure  arrays  by  SBLDO  or  SBLD01. 

NBLD01  (I , J , V,MD,  FD , RY , CY , A,  AN,  NDR.NDC) 

Insert  element  (I,J)  into  the  sparsity  structure  arrays  of  a 
nonsymmetric  matrix  and  a corresponding  value  V into  the  coeffi- 
cient array. 

NBLD1  (I , J,V,MD,FD,A,  AN) 

Associate  a value  V to  element  (I,J)  of  a nonsymmetric  matrix. 
The  V-values  must  be  in  the  same  order  in  which  the  (I , J)  -values 
were  read- in  during  prior  construction  of  the  corresponding 
sparsity  structure  arrays  by  NBLD01. 


DEC  - Routines 
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SDECO(MD,RY,CY,ND, IP, IE,  IH) 

Generate  a record  of  the  decomposition  of  a symmetric  matrix 
on  the  basis  of  the  given  sparsity  structure. 
SDEC01(MD,FD,RY,CY,A,AN,ND,IP,IE,IH) 

Decompose  a given  symmetric  matrix  and  optionally  (MD(3)^0) 
generate  a record  of  the  decomposition. 

SDEC1(MD,FD,A,AN,IE) 

Decompose  a given  symmetric  matrix  using  a previously  generated 
decomposition  record. 

NDEC01(MD,FD,RY>CYfA,AN,NDR,NDC,IPR,IPC,IE,IHfNGl,Nu2) 

Decompose  a given  nonsymmetric  matrix  and  optionally  (MD(3)^0) 
generate  a decomposition  record. 

NDEC1(MD,FD,A,AN,IE,IH) 

Decompose  a given  nonsymmetric  matrix  using  a previously  generated 
decomposition  record. 

Pivot  Routine 

PVT01  (I , N,  IX , KR,  KC , F , RY , CY , A , IPR,  IPC , NDR , NDC , IE , IH) 

Select  the  next  pivot  by  the  minimal  degree  algorithm  during  the 
decomposition  of  a nonsymmetric  matrix.  The  routine  is  used  by 


NDEC01 . 
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SLV  - Routines 

These  routines  use  the  decomposed  matrix  in  secondary  storage. 

The  right  side  of  the  system  is  given  in  the  form  of  the  input  array 
X which  in  turn  may  be  the  same  as  the  output  array  Y of  the  solu- 
tion. The  routines  may  be  called  repeatedly  for  different  right  sides. 
SSLV(MD,X,Y) 

Return  the  solution  Y of  the  symmetric  system  with  right 
side  in  X. 

NSLV(MD,X,Y,AN) 

Return  the  solution  Y of  the  nonsymmetric  system  with  right 
side  in  X. 

I/O  - Routines 

The  I/O  routines  for  communication  with  the  random  access  storage 
device  are  not  in  basic  FORTRAN.  They  should  be  modified  to  suit  a user's 
machine  configuration.  The  routines  have  the  following  entries: 

For  I/O  of  decomposition  record  (array  of  positive  integers) 

DWI  - Initialize  for  writing. 

DW(K)  - Write  K as  next  entry  of  the  array. 

DWE  - Terminate  writing. 

DRI  - Initialize  for  reading. 

DR(K)  - Return  the  next  entry  of  the  array  in  K. 

For  I/O  of  symnetric  decomposed  matrix  (array  of  pairs) 

SVWI  - Initialize  for  writing. 

SVW(£,s)  - Write  (£,s)  as  next  entry  of  the  array,  ^-signed  integer, 


s-real. 
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SVWR  - Terminate  writing. 

SVRI  - Initialize  for  reading  (forward). 


SVRF(£,s)  - Return  the  next  entry  of  the  array  in  l and  s. 
SVRB(£,s)  - Return  the  previous  entry  of  the  array  in  l and  s. 
For  I/O  of  a nonsymnetric  decomposed  matrix  (two  arrays  of  pairs) 
NVWI  - Initialize  both  files  for  writing. 

NVWF(£,s)  - Write  (£,s)  as  next  entry  of  the  L-array  (^-signed 
integer,  s-real). 

NVWB(£,s)  - Write  (£,s)  as  next  entry  of  the  U-array. 

NVWE  - Terminate  writing  of  both  files. 

NVRI  - Initialize  for  reading,  file  L forward,  file  U 
backward . 

NVRF(«,s)  - Return  next  entry  of  L-array  in  (£,s). 

NVRB(£,s)  - Return  previous  entry  of  U-array  in  (£,s). 

The  following  Table  1 shows  the  usage  of  the  I/O  routines  by  the 
various  main  routines: 


Table  1 

I/O  Routine  Usage 

DNR  SVWR  NVWR 


Decomposition  Symmetric  Nonsytrmetric 
Record  Dec.  Matrix  Dec.  Matrix 

File  10  FUe  11  Files  12,  13 


Progrta 

Write  | Read 

Write  | Read 

Write  | Read 

SDECO 

X 

l 

| 

1 

SEBC01 

t . 

X - 

1 

seen 

1 X 

X | 

1 

SSLV 

1 

1 X 

.1 

roecoi 

t | 

1 

X | 

Ncea 

1 X 

1 

X | 

NSLV 

1 

1 X 

X • routine  used 

t * use  Is  optional,  depending  on  user's  request 
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4.2  Arguments : The  arguments  in  the  various  calling  sequences  either 
reference  single  values  or  data  arrays.  For  simplicity  the  single  varia- 
bles are  collected  in  two  arrays,  an  integer  array 

MD(I),  I = 1,2,. ..,8 

and  a real  array 

FD(I) , I = 1,2,. ..,7  . 

The  first  three  values  of  MD  and  the  first  two  of  FD  are  to  be  supplied 
by  the  user;  the  others  represent  output  of  various  other  routines.  Care 
should  be  taken  that  these  latter  values  are  not  modified  whenever  they 
are  still  to  be  used  as  input  by  other  routines. 

The  use  of  the  various  arguments  by  the  routines  in  the  package  is 
summarized  in  Tables  2 and  3 below. 

MD  - Array 

MD(1)  = N The  dimension  of  the  matrix;  to  be  supplied  by  the  user. 

MD(2)  = MX  The  lengths  of  the  arrays  RY,  CY  and  A.  If  the  decomposi- 
tion of  the  matrix  requires  more  internal  storage,  that 
is,  if  > MX,  then  the  error  indicator  MD(4)  is  set  to 
one  and  the  process  is  terminated  with  a return  to  the 
user's  main  program. 

MD(3) 


If  this  indicator  is  zero,  the  DEC01  program  does  not 
produce  a decomposition  record;  for  any  nonzero  value 
of  MD(3)  such  a record  is  generated. 
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MD(4)  Error  indicator  set  as  follows: 

= 0 : no  error. 

= 1 : storage  overflow;  MX  is  too  small  for  decomposition. 

= 2 : on  the  basis  of  the  sparsity  structure  (independent 

of  the  values  of  the  elements)  the  matrix  is  singular. 
= 3 : the  matrix  is  declared  numerically  singular. 

MD(5)  In  the  nonsymmetric  case  equal  to  MO  + N where  MO  is  the 
number  of  nonzero  elements  in  the  matrix;  in  the  symmetric 
case  equal  to  MO,  the  number  of  nonzero  elements  on  the 
diagonal  and  in  the  upper  (or  lower)  triangle  of  the  matrix. 
MD(6)  The  length  of  the  actually  utilized  portion  of  the  arrays 
RY,  CY  or  A,  equal  to  ML  + N or  Ml  in  the  nonsymmetric  or 
symmetric  case,  respectively. 

MD(7)  In  the  nonsymmetric  case  equal  to  M2  + N where  M2  is  the 

number  of  nonzero  elements  in  the  decomposed  matrix;  in  the 
symmetric  case  equal  to  the  nonzero  elements  on  the  diagonal 
and  in  the  lower  triangle  after  the  decomposition. 

MD(8)  Length  of  the  decomposition  record. 

FD  - Array 

FD(1)  A tolerance  EPS  supplied  by  the  user.  If  a pivot  value  in 
magnitude  is  less  than  EPS  the  matrix  is  considered  to  be 
numerically  singular  and  the  decomposition  is  terminated. 
FD(2)  Initial  input  by  the  user  containing  the  pivot  selection 
parameter  p. 

FD(3)  Largest  coefficient  value  in  magnitude  in  the  original 

matrix.  Initialized  by  the  IMT  routines  and  updated  by  the 
BLD  routines. 
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FD(4)  Largest  coefficient  value  in  magnitude  encountered  during 
decomposition,  calculated  by  the  DEC  routines. 

FD(5)  Natural  logarithm  of  the  absolute  value  of  the  determinant 
calculated  by  the  DEC  routines. 

FD(6)  The  sign  of  the  determinant  as  +1.0  or  -1.0  calculated  by 
the  DEC  routines. 

FD(7)  The  natural  logarithm  of  the  product  of  the  L2  norms 
of  the  row  vectors  of  the  original  matrix  calculated  by 
the  DEC  routines. 


Data  Arrays 


RY  (MX) , CY  (MX) 
A (MX) 

AN(N) 

ND(N) 

NDR(N)  ,NDC(N) 


IE(N),IH(N) 

IP(N) 

IPR(N) , I PC (N) 
NG1(N),NG2(N) 


Integer  arrays  for  the  sparsity  structure.  Their 
length  MX  is  specified  by  MD(2). 

Real  array  for  the  values  of  the  nonzero  matrix 
elements. 

Real  array  of  length  N (see  MD(1))  used  to  collect 
row- vector  norms  of  the  matrix. 

For  symmetric  A. 

For  nonsynBetric  A.  Integer  arrays  of  length  N 
containing  the  number  of  nonzero  elements  by  row 
(or  column) . 

For  any  matrix. 

For  symmetric  matrices. 

For  nonsyrrmetric  matrices.  Temporary  arrays  of 
length  N. 


9 - 

Ibccept  for  the  last  seven  temporary  arrays  all  data  arrays  are 
initialized  in  the  appropriate  IKV  routines.  All  integer  arrays  are  used 
only  for  storing  nonnegative  integers.  Thus  for  particular  computers 
these  arrays  could  be  packed  together. 
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5.  Some  Experimental  Results 

Two  groups  of  computational  experiments  were  conducted  on  the  Univac 
1108  of  the  University  of  Maryland,  Computer  Science  Center.  They  correspond 
to  the  computational  experiments  reported  in  [1].  Since  the  basic  decomposi- 
tion procedure  is  the  same  as  in  [1  ],  the  overall  elapsed  time  for  execution 
of  the  decomposition  programs  (Tyj)  and  the  number  of  elements  after  decomposi- 
tion (Mp)  are  essentially  the  same  as  reported  there. 

The  new  results  presented  below  concern  the  maximal  in-core  storage 
requirement  (M.),  the  elapsed  time  for  execution  of  the  decomposition  routines 
(T^)  using  a previously  generated  decomposition  record,  and  the  elapsed  time 
for  execution  of  the  backsubstitution  routines  (T^)  when  the  decomposed 
matrices  are  residing  on  auxiliary  storage.  These  times  are  given  below 
relative  to  the  elapsed  time  (T^y)  for  the  execution  of  the  decomposition 
programs.  It  should  be  pointed  out  that  for  N larger  than  100  there  is 
less  than  a three  percent  difference  between  the  elapsed  time  for  the  genera- 
tion of  a decomposition  record  (SDECO)  and  that  for  a regular  decomposition 
with  or  without  retaining  the  record  (SDEC01,  NDEC01). 

The  first  group  of  experiments  involved  nonsymmetric  matrix  decomposi- 
tions. As  in  [1],  a special  program  was  used  to  generate  permuted  diagon- 
ally dominant  random  sparse  matrices  B = (b  „ ) . Given  the  dimension  N of 
B and  a number  Mq  of  nonzero  elements,  the  program  randomly  generates 
Mq  - N distinct  index  pairs  (i,j),  i ^ j,  1 s i,j  2 N,  and  the  corresponding 
matrix  elements  b^j . Then  each  diagonal  element  is  obtained  by  adding  a 
random  positive  number  to  the  sum  of  the  moduli  of  the  off-diagonal  elements 
in  its  row.  Finally,  the  rows  and  columns  are  independently  and  randomly 
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permuted.  Table  4 contains  the  results  obtained  when  the  decomposition 
programs  are  applied  to  these  random  matrices.  The  pivot  selection 
parameter  u = 0.125  was  used. 


Table  4 


N 

Results  for  Nonsymmetric  Random 

Matrices 

M 

0 

Mi 

YTw 

tsq/tlh 

50 

200 

200 

375 

.255 

.068 

50 

300 

300 

632 

.218 

.041 

100 

400 

400 

811 

.198 

.044 

100 

600 

1,114 

2,026 

.140 

.016 

200 

800 

913 

2,312 

.138 

.020 

200 

1,200 

3,991 

6,439 

.074 

.005 

300 

1,200 

1,924 

4,224 

.121 

.012 

The  second  group  of  experiments  involved  the  decomposition  of 
symmetric  matrices  obtained  by  the  discretizing  of  Dirichlet's  problem 
for  Laplace's  equation  on  the  unit  square  with  the  five-point  and  nine- 
point  formulas.  In  all  cases  a uniform  mesh  was  used.  The  coefficients 
of  the  resulting  matrices  are  well  known  and  are  not  repeated  here. 

Table  5 contains  the  results  obtained  with  the  symmetric  decomposition 


•programs. 


- 23  - 


5-point  Formula 

Table  5 

Symmetric  Matrix  Decomposition 

M0  M1  M2  VTLU 

TS0/TL! 

N = 81 

225 

225 

469 

.47 

.18 

N = 289 

833 

883 

2,413 

.22 

.06 

N = 484 

1,408 

1,699 

4,733 

.17 

.04 

N = 625 

1,825 

2,328 

6,566 

.16 

.03 

9-point  Formula 

N = 81 

353 

353 

715 

.35 

.10 

N = 289 

1,345 

1,673 

4,296 

.15 

.03 

N = 484 

2,290 

2,928 

8,054 

.13 

.02 

N = 625 

2,977 

4,047 

11,800 

.10 

.01 

6.  Program  Listings 
6. 1 Main  Package: 
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Ool 

0U2 

003 

004 
00b 
006 

007 

008 

009 

010 

011 

012 
01b 

014 

015 


SUBROUTINE  SINTO(MDfRYfCYfND) 
DIMENSION  MD(1 ) »RY(1) »CY(1) »ND(1 ) 
INTEGER  RYfCY 

C ***************************************** 

c * initialize  symmetric  structure  arrays  * 

C ***************************************** 

c 

N = MD ( 1 ) 

MD(b)  = N 
UO  10  I = 1fN 
RY(I)  = I 
CY(i)  = I 
lu  NO ( 1 ) = 1 
RETURN 
END 


001 

0U2 

003 

004 

005 

006 

007 

008 


SUbROUTlNE  SINT01(MDfFDfRYfCYfAfAN»ND) 

DIMENSION  MD(1)  ,FD(1)  |RY(1)'CY(1>  » A ( 1 ) , A,m  ( 1 ) , ND  ( 1 ) 
INTEGER  RYfCY 

C * INITIALIZE  SYMMETRIC  STRUCTURE  AND  COEFFICIENT  ARRAYS  ♦ 

C ********************************************************* 

N = MD(1) 


0U9 

010 

011 

012 

013 

014 

015 

016 
017 


MD(5)  = N 
FU(3)  = 0. 
DO  10  1 = 1 f N 
AN ( I ) = 0. 

R Y ( I ) = 1 
LY ( I ) - I 
10  NO ( 1 ) = 1 
RETURN 
END 


UOl 

002 

003 

004 

005 

006 

007 

008 

009 

010 

011 

012 

013 

014 

015 


SUBROUTINE  SINTI(MDfFDfAN) 

DIMENSION  MD(1)fFD(1)fAN(1) 

C ****************************************** 
C * INITIALIZE  SYMMETRIC  COEFFICIENT  ARRAY  * 
C ****************************************** 

C 

N = MD ( 1 ) 

MD(5)  = N 
FD  ( 3 ) = 0. 

MO  = MD(5> 

DO  10  1=1 rN 
10  AN ( 1 ) = 0 • 

RETURN 

END 

c 
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U 01 
UU2 
UU3 
U04 
UU5 
Uuo 
UU  7 
uua 
0 09 
010 
Oil 

SIS 

014 

01b 

01b 

017 

Old 


SUBROUTINE  NINTUl <MD»FD»RY*CY, AN,NDR,NDC) 

DIMENSION  MD{  1 ) »FD< 1 ) »RY ( 1 > »CY( 1 ) t NDR  ( 1 ) ,NUt(  1 ) * AN<  1 ) 
INTEGER  R Y * C Y 

C * INITIALIZE  NONSYMMETRIC  STRUCTURE  AND  COEFFICIENT  ARRAYS  * 
C 

N - MD ( 1 ) 

MU ( b ) = N 
F U ( 3 ) - 0. 

00  10  1=1 # N 

SyH!  i !• 

nSr{{  ) “=*0 

lu  NDCll)  = 0 
RETURN 
END 


I 

I 


JOl 

002 

003 

004 
OUb 
00b 
Oo  7 
Oud 

009 

010 

011 

012 
013 


C 

c 

c 

c 


10 


SUBROUTINE  NINT1(MD»FD»AN) 

0 1 MENS  I ON  MD ( 1 ) #FD ( 1 ) » AN ( 1 ) 

♦ initialize  nonsymmetric  coefficient  array  * 


N = MU ( 1 ) 

MU ( b ) = N 
Fl)(3)  = U. 
UU  10  I = 1 » N 
AN ( 1 ) = 0. 

RETURN 
LNU 


| 


I 


u U 1 
002 

003 

004 
0 05 
00b 
007 

ooa 

009 

010 

011 

012 

013 

014 

015 
0 10 

017 

018 

019 

020 
021 
022 

023 

024 


C 

C 

C 

C 


10 


SubROUT  1 nE  SbLDU ( I * U t MD t R Y > C Y » ND ) 
U I MENS  ION  MD ( 1 ) rRY ( 1 ) »CY ( 1 ) » NO ( 1 ) 
INTEGER  KYiCY 

* HU1LU  SYMMETRIC  STRUCTURE  ARRAY  * 


IF  (I.EU.J)  RETURN 
MU  ( b ) = MU  ( b ) + 1 
My  = MQ(5) 

ND ( I ) = ND ( I ) + 1 

ND(J)  = NL) ( J ) + 1 
IF  ( I .GT. J)  GO  TO  10 
KY (MO ) = RY ( I ) 

CY(MO)  = CY(J) 

R Y ( I ) = MO 
CY(J>  = MU 
RETURN 

RY  ( MO  ) = R Y ( J ) 

CY (MU ) = CY ( I ) 

ry<u>  = mu 

CY ( i ) = MU 

RETURN 

END 
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uui 

0U2 

U03 


UU4 

OUb 

006 

0U7 

000 

U09 

010 

011 

012 
J 13 
014 
01b 
016 


017 

oia 

019 

020 


SUBROUTINE  SBLD01 ( I »J#V#MD>FD,RY»CY, A»AN,ND) 
DIMENSION  MD( 1 ) »F0<1) »RY(1 ) »CY(1 ) * A(1 ) » AN(1 ) »ND(1 » 
INTEGER  R Y » C Y 

C **************************************************** 

C * BUILD  SYMMETRIC  STRUCTURE  AND  COEFFICIENT  ARRAYS  * 

C **************************************************** 

c 

b = V** 2 

FD ( 3 ) = AMAX1 ( FD  l 3 ) i ABS » V ) ) 

AN  ( 1 ) = AN  ( I ) +S 
IF  (I.EQ.J)  GO  TO  20 
MO(b)  = MU ( b ) 4 1 
MO  = MD(5) 

A(MO)  = V 
AN ( J ) = AN ( J ) +S 
ND( I ) = ND ( I ) +1 
ND(J)  = NU ( J } + 1 
IF  (I.GT.J)  GO  TO  10 
KY (MO ) = K Y ( I ) 

CY (MO)  = CY ( J) 


021 

U22 

023 

024  10 
02b 

02b 

027 

028 

U29  20 

030 

031 


KY ( I ) = MO 
CY(J>  = MO 
RETURN 

RY(MO)  = R Y ( J ) 
CY(MO)  = CY  ( I ) 
RY (U)  = MU 
CY ( 1 ) = MO 
RETURN 
A ( I ) = V 
RETURN 
END 


001 

002 

003 

004 
00b 
006 
007 
000 

009 

010 

011 

012 

013 

014 

015 

016 

017 

018 


SUBROUTINE  SdLDl( I * Jr  V , MU , FD» A # AN ) 
UIMENSION  MDC1)  »FD(1)»A(1)»ANU) 

C ************************************* 

C * BUILD  SYMMETRIC  COEFFICIENT  ARRAY  * 

C ************************************* 

c 

b = v**2 

FD(3)  = AMAX1 (FD<3) i ABS(V) ) 

AN ( I ) = AN ( I ) +S 
IF  (I.EQ.J)  GO  TO  lU 
MD(b)  = MD(b)+l 
MO  = MD(b) 

A ( Mo ) = V 

AN ( J ) = AN ( J ) +S 

RETURN 

10  A ( I ) = V 

RETURN 

END 
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uol 
002 
003 
J 04 
OOb 
00b 
0 07 
008 

009 

010 

011 

012 

013 

014 
01b 
016 

017 

018 

019 

020 


C 

C 

C 

C 


SUBROUTINE  NBLD01  ( I#J»V»MD*FD,RY,CY# A» AN,NDR#NDC) 

DIMENSION  MD(l)»rD(l)*RY(l)#CY(l)#A(l ) t NOR ( 1 ) * NDC  ( 1 ) » AN  ( 1 ) 
INTEGER  RY  » C Y 

* BUILD  NONSYIMMETRIC  STRUCTURE  AND  COEFFICIENT  ARRAYS  * 


MD(b)  = MU ( 5 ) + 1 
MO  = MD(b) 

KY(MO)  = R Y ( I ) 

OY(MO)  = CY(U) 

K Y ( I ) = MO 
Cy(J)  = Mo 
NUR(I)  = NDR ( I ) +1 
NDC(J)  = NDC ( J ) + 1 
A ( MU ) = V 

AN  ( 1 ) = AN(I)+V**2 

Eu<3)  = AMAX1(FD(3)»ABS(\/)) 

Kfc Turn 

tNO 


001 

002 
0 03 


SUBROUTINE  NBLOK  I > J»  V » M() , FD»  A * AN ) 
U I MENS  I UN  MU(  1 ) #F0(1 ) * A(  1 ) » AN(  1 ) 


004 
0 0b 
006 

007 

008 


009 

010 

011 

012 
013 


C 

C 

L 


* BUlLU  NONSYMMETRIC  COEFFICIENT  ARRAY  ♦ 


AN ( 1 ) z AN ( I ) +V**2 
MD(b)  = MU ( b ) +1 
Mo  = M U ( b ) 

A'Mu)  = V 

FD ( 3 ) = AMAX1 ( FD ( 3 ) » ABS ( V ) ? 

RETURN 

ENU 


oul 


002 

003 

004 


00b 

006 

0U7 

008 

009 

010 

011 

012 

013 

014 


01b 
U 16 
017 


SUBROUTINE  SOECO(MDiRY»CY»ND» IPr IE» IH) 

DIMENSION  MD(1)  »RY(1)#CY(1)»ND(1),IP(1),IEU)#IH(1) 
INTEGER  RY ♦ C Y 

C ♦ GENERATE  SYMMETRIC  DECOMPOSITION  RECORD  * 

C **************  **********************  ******* 


MM  = 0 
N = MO ( 1 ) 

MU (4)  = 0 
MQ(b)  = MU ( b ) 

MO  ( 7 ) = MD(b) 

MO ( 8 ) = 0 

C INITIALIZE  AUXILIARY  FILE  FOR  WRITING 
CALL  Owl 
DO  10  I - 1 » N 

lu  IPU)  = I 


r 
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I 


U 18 

019 

U20 

021 

022 

023 

024 
02b 
026 

027 

028 

029 

030 

031 
0 32 

033 

034 
03b 

036 

037 


C 

C 

C 


c 

c 

c 


20 


LOOP  ON  PIVOTING 

UO  2b0  I=1#N 
K = 0 

IF  (I.EQ.N)  GO  TO  30 

SELECT  PIVOT  BY  MINIMAL  DEGREE 

NDX  = N+l 
UO  20  J-I » N 
IX  = IP(J) 

IK  (i\lD(  IX)  .GE.NDX)  GO  TO  20 
NDX  = NU(1X) 

1Y  = J 
CONTINUE 

IF  (I.EQ.1Y)  GO  TO  30 
J = IP( 1Y) 

IP ( 1Y  ) = IP ( I ) 

IP ( I ) = J 


030 

C 

039 

C 

COLLECT  THE  ROW  AND  COLUMN 

040 

C 

ALSO  DELETE  THEM  FROM  TH 

041 

C 

042 

30 

IX  = IP l I ) 

k 

043 

IF  (I.EQ.N)  GO  TO  HO 

044 

I Y 1 = 0 

04b 

IY  = IX 

04b 

40 

IY  = RY(IY) 

047 

IF  (IY.EO.IX)  GO  TO  70 

I 

048 

IZ  = IY 

049 

bu 

IZ  = CY(IZ) 

ObO 

Obl 

Ob2 

0b3 

Ob4 

0b5 

Ob6 

0b7 

Ob8 

Ob9 

ObO 

Obl 

002 
0b3 
064 
Uob 
Obb 
Ub7 
Ob8 
Ob9 
0 70 
071 
)72 

073 

074 
07b 
0 7b 
U 77 

078 

079 
000 
081 
0b2 

003 
084 
OOb 
006 

007 

008 


IF  (IZ.GT.N)  GO  TO  60 
K = K + 1 
IE(K)  = IY 
I H ( K ) = IZ 
ND(IZ)  = ND(IZ)-1 
bU  IF  ICY ( 1Z ) .NE. IY)  GO  TO  50 
CY(1Z)  = CY(lY) 

CY!lr)  = mm 
MM  = IY 
oO  TO  40 

7 0 IY  = CY(IY) 

IF  (IY.EU.IX)  GO  TO  100 
iZ  = IY 
I Y 1 = IY 

00  IZ  = KY(IZ) 

IF  (IZ. GT,N)  GO  TO  90 
K = K+l 
IE ( K ) = IY 
IH(K)  = IZ 
ND(1Z)  = ND  ( IZ  ) - 1 
90  IF  (KY(IZ).NE.lY)  GO  TO  80 

K Y ( IZ)  = RY( IY) 

GO  TO  70 

100  IF  (1Y1.EQ.0)  GO  TO  110 
CY(IYl)  = MM 
MM  = CY(1X) 

C 

C MODIFICATION  OF  THE  ROw  ELEMENTS 
C 

C WRITE  OUT  PIVOT 

110  CALL  Dw(IX) 

CALL  UW(K) 

MD(O)  = MU ( 8 ) ♦ 1 +K**2 
IF  (K.EO.O)  GO  TO  2b0 
DO  lib  J=1*K 
CALL  DW(IE( J)  ) 
lib  CALL  DwUH(  J)  ) 

IF  (K.EQ.l)  GO  TO  250 

R 1 = K-l 


'r 


W 
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069 

USJO 

091 


112 

lii 

114 

115 
no 
117 
lid 

119 

120 
HI 
122 

123 

124 
12b 

1 2b 
127 
I2d 
129 
liQ 
lil 

132 

133 

134 
lib 
lib 
137 
13Q 

139 

140 

141 

142 

143 

144 
14b 

146 

147 
14d 
149 
lbO 
lbl 


c 

c 

c 


LOOP  FOR  THE  CROSS-POINT  ELEMENTS 


0 92 

UO 

240  J- 1 * K 1 

093 

Jl 

— J + 1 

U94 

u 

= IFI  ( J ) 

09b 

Do 

230  0J=J1»K 

096 

02 

= IH(JJ) 

097 

il 

= M INO ( I Z » JZ ) 

09ti 

1 2 

= MAX  U ( IZ» JZ) 

099 

LI 

= RY(I1) 

100 

L2 

= C Y < 12 ) 

101 

120 

IF 

( ( L 1 , EU , I 1 ) .OR 

102 

IF 

(Ll.EQ.L2)  GO 

loi 

IF 

(L1.GT.L2)  GO 

1 04 

L2 

= C Y ( L2 ) 

10b 

OU 

10  120 

106 

1 3U 

LI 

= RY(L1) 

107 

GO 

10  120 

1 06 

c 

INSERTION  OF  A NEW  N 

109 

140 

ND(ll)  = NU ( Ill+l 

no 

NUII2)  = NU ( 1 2 ) + 1 

ill 

MD(7)  = MU(7)+1 

1 b 0 


IF  (MM.NE.O)  GO  TO  170 
USE  NEw  ST  OR AGE 
MU ( b ) = MO ( 6 ) ♦ 1 
IF  <MD(6) . LE ♦ MD ( 2 ) ) 

MD ( 4 ) = 0 
RETURN 
LI  = MD(b) 

RY(L1)  = RY(ll) 

CY(L1)  = C Y < 12) 

RYdl)  = LI 
CY ( 12 ) = LI 


GO  TO  160 


C USE  AVA 
170 


oO  TO  220 
ILAbLE 


STORAGE 


LI  = MM 
MM  = CY(MM) 

L3  = II 
L2  = 12 

IdO  IF  (RY(L3) .LT.L1)  GO  TO  190 
L3  = RY(L3) 

GO  TO  160 

190  RY(L1)  = RY(L3) 

RY(L3)  = LI 

2U0  IF  (CY(L2) .LT.Li)  GO  TO  210 
L2  = CYIL2) 

GO  TO  200 

210  CY(L1)  = CY(L2) 

Cy(L2>  = LI 

C WRITE  OUT  CROSS-POINT  ELEMENT 
22 U CALL  DW(L1) 

230  CONTINUE 

240  CONI INUE 

C END  OF  MODIFICATION  LOOP 

2b0  CONIINUE 

C 

C END  OF  PIVOTING  LOOP 
C 

CALL  DwE 
RETURN 
C 


ENO 


30  - 


ool 
U U 2 
0 03 
OU4 
005 
U 06 

007 

008 

009 

010 
Ull 
012 
U 1 3 
014 
01b 
016 

017 

018 

019 

020 
021 


SUBHOUT  1NE  SDEC01(MU,FD»«Y#CY»A,AN»ND,  IP,  IE,  1 H ) 

DIMENSION  M0(1 ) ,FD(1 ) »RY  < 1 ) »CY(1 ) »NU(1) , IP( 1 ) i IE(1)  ,IH(D 
DIMENSION  A(1),AN(1) 

INTEGER  R Y # C Y 

C ********************************************* 

C ♦ DECOMPOSE  SYMMETRIC  MATRIX  AND  OPTIONALLY  * 

C ♦ GENERATE  DECOMPOSITION  RECORD  * 

C ********************************************* 

c 

N = MD ( 1 ) 

MD(4)  = 0 
FD  ( b ) = 0 . 

FD  ( 6 ) = 1. 

FD(7)  = 0. 

FD(4)  = 0. 

MM  = 0 

MD(b)  = MD(b> 

MD(7)  = MD(S) 

MO  (8)  = 0 

C INITIALIZE  AUXILIARY  FILE  FOR  WRITING 
IF  1 MD ( 3 ) .NE.O)  CALL  DW  I 


022 

023 

024 
02b 
026 

027 

028 

029 

030 

031 
0 32 

033 

034 
03b 
0 36 

037 

038 

039 

040 

041 

042 

043 

044 

045 

046 

047 

048 

049 
0 50 
Obi 
0 52 
053 
0 54 

055 

056 

057 

058 

059 

060 
Obi 
062 


DO  10  I=1»N 

FD(7)  = FD ( 7 ) + ALOG ( AN  { I ) ) 

1U  I P ( I ) = I 

FD ( 7 ) = 0 . 5*FD ( 7 ) 

CALL  SVwi 

C 

C LOOP  ON  PIVOTING 
C 

DO  250  I - 1 » N 
K = 0 

IF  (I.EQ.N)  GO  TO  30 
C 

C StLECT  PIVOT  BY  MINIMAL  DEGREE 
C 

NDX  = N+l 
AX  = 0. 

DO  15  J - 1 » N 
IX  = 1P(J) 

15  AX  = AMAX1(AX*ABS(A(IX) ) ) 

AX  = AX*FD(2) 

DO  20  J-I » N 
IX  = 1P(J) 

IF  (ABS(A( IX) ) .LT.AX)  GO  TO  20 
IF  (ND(iX) .GE.NDX)  GO  TO  20 
NDX  = ND ( I X ) 

1Y  = J 

2U  CONTINUE 

IF  (I.EQ.lY)  GO  TO  30 
J = IP ( 1 Y ) 

IP ( I Y ) = IP ( I ) 

IP(1)  = J 
C 

C COLLECT  THE  ROW  AND  COLUMN  OF  ThE  PIVOT 

STORAGE 


TO  300 


IY  = IX 


C 

C 

3u 


ALSO  DELETE  THEM  FROM  THE 

IX  = IP ( I ) 
b = A ( I X ) 

IF  <AB$(S) .LT.FDI 1 ) ) GO 
IF  (I.EQ.N)  GO  TO  1 1 0 
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I 


l 


4|J  1Y  = RY  l It  ) 

IF  (IY.LU.IX)  GO  TO  70 

u - ir 

bu  n = c Y ( IZ) 

IF  (IZ.GT.N)  GO  TO  bO 
K = K + l 
I£(K)  = IY 
1H ( K ) = ll 
AN ( K ) = A ( I Y ) 

A ( I Y ) = U . 

ND( IZ)  = NDl IZ)-1 
UU  IF  ICY(IZ).NE.IY)  GO  TO  bO 

CYUZ)  = CY(IY) 

CY(IY)  = MM 
MM  = IY 
GO  10  4U 

7 0 IY  = CY(IY) 

IF  (1Y.L0.IX)  GO  TO  100 
ll  - IY 
1 Y 1 = IY 

au  ll  = RYUZ) 

IF  (IZ.GT.N)  GO  TO  90 
K = K + l 
IL(K)  = IY 
i H ( K ) = ll 
AN ( K ) = A < I Y ) 

A ( I Y ) = U . 
i\ID(1Z)  = NQ(iZ)-l 
9u  IF  (RY(  1Z)  ,n£,  IY)  GO  TO  BO 

KY(1Z)  = KY(lY) 

00  10  7U 

iuu  IF  (lYl.tU.U)  GO  TO  110 

CYUYl)  = MM 
MM  r C Y ( I X ) 

C 

C MODIFICATION  OF  lHt  ROw  ELEMENTS 
C 

C fcRllt  OUT  PIVOT 

11U  IF  ( MD ( 3 ) . NE  . 0 ) CALL  DW(IX) 

IF  ( MD { 3 ) , NL , 0 ) CALL  DW ( K ) 

FD ( 4 ) = AMAX1 (FD(4) » ABS(S)  ) 

F D ( b ) = FU(b)+ALOG( ABS(S)  ) 

IF  (S.LT.O. ) FD ( 6 ) = -FD ( 6 ) 

CALL  SVW(-IXiS) 

MU ( a ) = MU( 8 ) +1+K**2 
IF  (K.EO.O)  GO  TO  2b0 
UO  115  J=lrK 
AN( J)  = AN(J)/S 
CALL  SVW( IH( J) » AN ( J ) ) 

F0(4)  = AtoAXl ( FD ( 4 ) » ABS ( AN ( J) ) ) 
IF  ( ML)  ( 3 ) • NE  , 0 ) CALL  DW  ( I E ( J ) ) 
lib  IF  ( MU ( 3 ) • NE . 0 ) CALL  UW<IH(U)) 

K 1 - K— 1 
C 

C LOOP  Fur  the  cross-point  elements 

c 


DO  240  0=1, K 


- J+ 1 
= IFM  J ) 
AN  ( J ) 


A ( IZ ) = A( IZ)-S*Z**2 


FD ( 4 ) = AMAX1(FD(4) »ABS(A(IZ) ) ) 


IF  (J.EO.K)  GO  TO  240 
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12o 

127 

126 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 
14  0 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 
io9 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 

187 

188 

189 

190 


DO  230  JJ=J1»K 
JZ  = IH(JJ) 

11  = MINO(IZ'JZ) 

12  = MAX0(IZ»JZ) 

LI  = RY(I1) 

120  T(Li!Iq!xi ) ,0R, (L2.E0.I2)  ) GO  TO  140 

IF  (Ll,tU,L2)  GO  TO  220 
IF  (L1.GT.L2)  GO  TO  130 
L2  = CYIL2) 

GO  TO  120 
130  LI  = RY(L1) 

GO  TO  120 

C INSERTION  OF  A NEW  NON-ZERO  ELEMENT 
1 4 0 ND (11)  = ND ( 1 1 ) +1 
NO ( 12 ) = ND ( 12 ) +1 
MD(7)  = M0{7)+1 
IF  (MM.NE.O)  GO  TO  170 
C USE  NEW  STORAGE 

MU (6)  = MU ( 6 ) + 1 
IF  ( MD (6 ) ,GT . MD (2 ) ) GO  TO  310 
lbO  LI  = MD ( 6 ) 

A ( L 1 ) = 0. 

KY (LI ) = RY ( II ) 

CY(L1)  = C Y ( 12) 

RY  U1  ) = LI 
CY(  12)  = LI 
GO  TO  220 

C USE  AVAILABLE  STORAGE 
170  LI  = MM 

MM  = CY(MM) 

L3  = 11 
L2  = 12 

IF  (RY(L3).LT.L1)  GO  TO  190 
L3  = RYIL3) 

GO  10  180 
RY(L1)  = RY(L3) 

RY(L3)  = LI 

IF  (CY(L2) .LT.L1)  GO  TO  210 
L2  = C Y ( L2 ) 

GO  TO  200 
CY(Ll)  = LY(L2) 

CY(L2)  = Li 
C WRITE  OUT  CROSS-POINT  ELEMENT 
220  IF  ( MD ( 3 ) , NE . 0 ) CALL  DW(L1) 

A(L1)  = A ( LI ) -AN ( J)*AN( JJ)*S 
230  Fu (4 ) = AMAX1(FD(4)pABS(A(L1))) 

240  CONTINUE 

C END  OF  MODIFICATION  LOOP 
250  CALL  SVW(-IXfS) 

C 

END  OF  PIVOTING  LOOP 


leO 


190 

200 


210 


C 

C 


CALL  SVWE 

IF  ( MD ( 3 ) , NE . 0 ) CALL  DwE 
RETURN 


C 

C SINGULAR  MATRIX 
C 

300 


310 

C 


MD ( 4 ) = 1 
RETURN 
MD(4)  = 3 
RETURN 

END 


/ ^ 


-ir~ 


T ‘ 
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001 

002 

U03 

004 

005 
OOfa 
007 
0U8 

009 

010 

011 

012 

013 

014 
01b 
016 

017 

018 

019 

020 
021 
022 

023 

024 
02b 
026 

027 

028 

029 

030 

031 
0 32 

033 

034 

035 

036 

037 

038 
0 39 
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04b 
04b 

047 

048 

049 
ObO 


SUBROUTINE  SUECl<MD,FD,A,AN,IE) 
DIMENSION  MD(l),FD<l)*A(l)fAN(l),lE(l) 


C 

N = MD ( 1 ) 

MD ( 4 ) - 0 
FD ( b ) = 0. 

FD ( b ) = 1. 

FD ( 7 ) = 0. 

FD ( 4 ) = 0. 

DO  10  1=1 i N 

lu  Fd ( 7 ) = FD(7)+AL0G(AN< I ) ) 

FD(7)  = 0 . b*FO ( 7 ) 

C CLEAR  SlORAGt  TO  BE  FILLED 

IF  (MU(b) .LE.ND(b) > GO  TO  30 
i*l M = MD  ( b ) +1 
Ml  = MD(6) 

DD  20  I=MM,M1 

£uinit  ZE  P5r  read-in  and  WRITE-OUT 
3o  Call  dri 

CALL  SVWi 

C 

C LOOP  ON  THE  PIVOTS 

C 

DO  90  1 = 1 * N 
CALL  DR (IX) 

CALL  OR  l K ) 
b = A( IX) 

IF  (ABS(S)  .LT.FDU  ) ) GO  TO  100 
FD  < 4 ) = AMAX1 (FD(4) » ABS(S) ) 

F u ( b ) = FD(5)+ALOG(ABS(S)  ) 

IF  (S.LT.O.)  FD ( 6 ) = -FD(6) 

IF  (ts.Eu.O)  GO  TO  7U 
C 

C COLLECT  THE  ROW  OF  THE  PIVOT 
C 

DO  40  J=1*K 
CALL  OR (LI) 

CALL  DR ( IE ( J ) ) 

AN ( J ) = A(L1 ) 

A ( L 1 ) = 0. 

4u  F U ( 4 ) = AMAX1(FD(4) »ABS(AN<J)  ) ) 

C 

C MODIFICATION  OF  THE  ROW  ELEMENTS 
C 

DO  bO  J=1»K 
1Z  = IE ( J) 


Obi 
0b2 
0b3 
0b4 
055 
0 56 

057 

058 

059 
ObO 
061 
062 
063 
0b4 
065 


Z = AN ( J ) 

AN(D)  = AM  J ) /S 
A ( I Z ) = A( IZ)-Z*AN( J) 

FD ( 4 ) = AMAX1 (FD ( 4 ) r ABS ( AN ( J) ) ) 

FD ( 4 ) = AMAX1(FD(4) »ABS(A(IZ) ) ) 

IF  (J.EO.K)  GO  TO  60 
J1  = J+l 
C 

C MODIFICATION  OF  THE  CROSS-POINT  ELEMENTS 
C 

DO  50  JJ=J1»K 
CALL  DR(L1) 

A ( LI  ) = A<L1)-AM  JUAM  JJ) 

50  FD  ( 4 ) = ANiAXl  (FD(4)  r ABS(A(L1)  ) ) 

bO  CONI INUt 
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066 

Ub7 

066 

069 

070 

071 

072 

073 

074 
0 7b 
07b 
0/7 

078 

079 
060 
061 
062 
083 
064 
085 


C WRITE  OUT  PIVOT  AMD  ITS  ROW 
C 

70  CALL  SV*(-IX,S) 

IF  (K.EO.O)  GO  TO  90 
UO  60  J=1»K 

60  CALL  SVW(  IE( J) »AN( J) ) 

90  CALL  S VW l - 1 X » S ) 

C 

C END  OF  PIVOTING  LOOP 
C 

CALL  SVWt 
RETURN 
C 

C SINGULAR  MATRIX 

C 

100  MO ( 4 ) = 3 

RETURN 

END 


001 

002 

003 

004 

005 

006 

007 

008 

009 

010 

011 

012 

013 

014 

015 

016 

017 

018 

019 

020 
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022 
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044 


SUBROUTINE  NDECU1(MD,FD,RY,CY, A , AN » NDR , NDC , IPR, IPC, IE, lH,NGl,NG2) 
DIMENSION  MO (1) > *0  < 1 ) , R Y ( 1 > » C Y ( 1 ) , A ( 1 ) , AN  < 1 ) » NDR ( 1 ) , NDC  ( 1 ) 
DIMENSION  IPR(1),IPC(1) »IE(1) , IH { 1 ) , NG1 ( 1 ) , NG2 ( 1 ) 

INTEGER  R Y » C Y 

C **DECOMPOSE*NONSYMMETRIC*MATR I X* AND*OPTIONALLY** 

c ********^*HJIE  ^EC0MPJSJTI0N  Record  ^ * 


10 


c 

c 

c 

c 

c 

c 


c 


MM  = 0 
N = MD ( 1 ) 

M D ( 4 ) = 0 
MD(b)  = iviU(5) 

MD ( 7 ) : MDI5I 
MD(8)  = 0 

IF  (MD(3) .NE.O)  CALL  DWl 
FD(5)  = 0. 

FD(b)  = 1. 

FD ( 7 ) = 0. 

F D ( 4 ) = 0. 

DO  10  1=1, N 

FD(7)  = FD ( 7 ) +ALQG ( AN ( I) ) 

IPR(I)  = i 
1PC(I)  = I 
FD ( 7 ) = 0.5*FD<7) 

Call  nvwi 

CALL  NVWb(-l,-l) 

LOOP  ON  PIVOTING 
DO  290  1=1, N 

PIVOT  SELECTION  BY  SEPARATE  PIVOTING  ROUTINE 

CALL  PVT01 ( I,M,IX,KR,KC,FD(2) ,RY ,CY , A, IPR, IPC,NDR 
IF  (KR.EQ.O)  GO  TO  310 
CHECK  PIVOT  VALUE  AND  UPDATE  DETERMINANT 
S = A ( IX ) 

IF  (ABS(S) .LE.FD(l) ) GO  TO  300 
FD ( 5 ) = FD(5)+AL0G(ABS(S) ) 

K = ( KR-fKC  ) /2 
K = KR+KC-2*K 
IF  (K.NE.O)  FD(b)  = -FD(6) 

FD ( 4 ) = AMAX1(FD(a) ,ABS(S) ) 


,NDC, IE, IH) 
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04b  c 

046  C COLLEC  1 THE.  ROW  AND  COLUMN  OF  THE  PIVOT 

047  C ALSO  FREE  THE  IK  STORAGE  LOCATIONS 
04b  C 

049  C THE  ROW 
UbO  M = U 

Obi  J = KR 

0b2  2u  J = RY(J) 

0b3  IF  (J.LE.N)  GO  TO  40 

0b4  IF  (J.EG.lX)  Go  TO  20 

Obb  M = Kl  + 1 

Obb  IE ( K 1 ) = J 

0b7  AN  ( M ) = A ( J ) 


0 b8 
0b9 
OdO 
Obi 
0t>2 


0b3 

Od4 

Obb 

Obb 

067 


ObB 


0b9 

070 

071 

072 

073 
0 74 
07b 


07b 

077 

u78 


079 

JbO 

Obi 

0b2 

0b3 

0M4 


Obb 

Obb 

0b7 

Obb 

0b9 

090 

091 

092 

093 

094 
U9b 

096 

097 

098 

099 

100 
101 
102 

103 

104 
10b 
i Ob 
107 


108 

109 

110 
111 


J 1 - J 

30  J1  = CY(J1) 

IF  (Jl.LE.N)  NG1 ( K 1 ) = J1 
IF  (Jl.LE.N)  NDC(Jl)  = NDC(J1)-1 
IF  ( C Y ( J1 ) . N£ . J ) GO  TO  30 
C Y ( J1 ) = C Y < J) 

L Y ( J ) = MM 
MM  = J 
GO  TO  20 
C THE  COLUMN 
4 U K2  - U 

J = KC 
K3  - 0 

bo  J = C Y ( J ) 

IF  (J.LE.N)  GO  TO  70 
K 3 - J 

IF  (J.EG.lX)  GO  TO  bU 
K2  = K2+1 
IH(K2)  = J 
A(K2)  = A < J ) /S 

Fu ( 4 ) = AMAX1(FD(4) *A8S(A(K2) ) > 

J1  = J 

oO  J1  = R Y ( J1 ) 

IF  (Jl.LE.N)  NG2(K2)  = J1 
IF  (Ji.Lt.N)  NDR(Jl)  = NDR(J1)-1 
IF  (RY  ( J1  ) , Nt , J ) GO  TO  bl) 

RY (ol ) = R Y ( J ) 
b 0 TO  50 

70  C Y ( K3 ) = NM 

MM  = CY (KC) 

C 

C WRITE  OUT  THE  PIVOT,  ITS  ROW  AND  COLUMN 
C AS  PARI  OF  (HE  DECOMPOSITION  RECORD 
C 

IF  ( MD ( 3 ) ,EQ,0 ) GO  TO  105 
CALL  Uw(IX) 

CALL  Dw(KC) 

CALL  DW(KR) 

CALL  DW(N2) 

ML ( b ) = MU(8)+(K1+1 )*(K2+1 ) 

IF  (K2.EO.O)  GO  TO  90 
DO  80  J=1»K2 
CALL  Dw ( IH( J) ) 
do  CALL  UW(NG2(J) ) 

9u  CALL  DW(Kl) 

IF  (Kl.EO.O)  GO  TO  105 
UO  100  J=1,K1 
CALL  Dw ( IE ( J) ) 

100  CALL  UW(NGKJ)  ) 

lob  IF  ( ( K l .tQ.O) .OR. (K2.EQ.0 ) ) GO  TO  230 
C 

C LOOP  TO  MODIFY  THE  INTERSECTING  ELEMENTS  BETWEEN  THE  ROW 
c and  column 

c 


I 


) 
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112 

113 

11*4 

lib 

116 

117 

116 

119 

120 
121 
122 

123 

124 
12b 

if  'i 

126 

129 

130 

131 

132 

133 

134 
13b 

136 

137 
136 

139 

140 

141 

142 

143 

144 
14b 

146 

147 
146 
149 
160 
lbl 
162 

163 

164 
166 
166 

167 

168 

169 
160 
161 
162 

163 

164 
16b 
166 

167 

168 
lb9 

170 

171 

172 

173 

174 
176 

176 

177 

178 


C 

110 


120 

c 

130 


DO  220  J=1,K2 
1Y  = NG2 ( J ) 

DO  220  K=1,K1 
K3  = RY  ( 1Y) 

OY  = NG1(K) 

K4  = C Y ( JYj. 

SEARCH  FOR  ELEMENT  IY»JY 
IF  (K3.EQ.K4)  GO  TO  210 
(K3.LT.K4)  GO  TO  120 
= R Y ( K3 ) 

(K3.LE.N)  GO  TO  130 
GO  TO  llu 
K4  = C Y ( K4  ) 

(K4.GT.N)  GO  TO  110 


IF 

K3 

IF 


IF 


IT  DOES  NOT  EXIST 
' ' = NDR 


NOR ( I Y j = NDR(1Y)+1 
NDC(JY)  = NDC  < JY ) + 1 
MO ( 7 ) = MD ( 7 ) +1 
IF  (MM.NE.O)  GO  TO  160 
C GET  NEW  LOCATION  FOR  THE  NEW  ELEMENT 
MD ( b ) = MD ( 6 ) ♦ 1 
IF  (MD(fa) ,LE,mO(2) ) GO  TO  150 

MD(4)  = 1 
RETURN 

160  K3  = MD(6) 

RY(K3)  = RY(IY) 

CY (K3>  = CY( JY) 

RY ( 1Y)  = K3 
CY ( JY)  = K3 
A(K3)  = U. 

GO  TO  210 

OLD  LOCATION  AVAILABLE  FOR  THE  NEW  ELEMENT 
K3  = MM 
A ( K3 ) = 0. 

MM  = CY(MM) 

K4  = IY 

IF  (HY  (K4) .LT.K3)  GO  TO  180 
K4  = RY(K4) 

GO  TO  170 
RY (K3)  = RY(K4) 

RY (K4  ) = K3 

IF  <CY(JY).LT.K3)  Go  TO  200 
JY  = CY(JY) 

GO  TO  190 
CY(K3)  = CY(JY) 

mod1^yY^lement 

IF  ( MD ( 3 ) , NE , 0 ) CALL  DW(K3) 

A(K3)  = A ( K3 ) "AN ( K ) * A ( J ) 

FD  ( 4 ) = AMAX1 (FD(4) »AB5(A(K3) ) ) 

220  CONTINUE 

C END  OF  MODIFICATION  LOOP 

C 

C WRITE  OUT  PIVOT  ROW  AND  COLUMN 
C 

CALL  NVWF(-KC»-KR) 

IF  (K2.EQ.0)  GO  TO  250 
DO  240  J=1»K2 
CALL  NVWF (NG2 ( J) # A ( J) ) 

IF  (Kl.EO.O)  GO  TO  270 
DO  260  J=1»K1 

260  CALL  NVWB(NG1(J)»AN(J)) 

270  CALL  NVWB(-KC»S) 

290  CONTINUE 
C END  OF  PIVOTING  LOOP 
C 

C END  FILES 


C 

lbO 


170 


160 

190 


200 

C 

210 


230 


240 

260 
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1 1 y 

L 

loo 

IF  ( MD { 3 ) , NE 

ldl 

CALL  NVWE 

lo2 

RETURN 

183 

C 

104 

C SINGULAR  MATRIX 

105 

C 

106 

300 

MD ( 4 ) = 3 

107 

RETURN 

188 

310 

MD ( 4 ) = 2 

ioy 

RETURN 

190 

C 

191 

END 

SUBROUTINE  PVT01( I*N, I X # KR * KC * F , R Y * C Y , A , i PH , IPC , NDR , NDC , IE, IH) 
DIMENSION  RY(l)»CY(l)#IPR(i)»IPC(l)#NDR(l)»NDC(l)»IE(l) , IH(] ) 
DIMENSION  A(l) 

INTEGER  RY,CY 

♦ MINIMAL  degree  pivoting  for  non-symmetric  matrix  * 

THE  routine*selects  pivot  by  minimal  degree*  it  is 

USED  BY  THE  ROUTINE  DECQl. 


IE  (l.NE.N)  GO  TO  30 
KR  = IPR(N) 

KC  = 1PC(N) 

IX  = RY  (KR) 

IF  (IX.NE.KR)  RETURN 
J KR  = 0 

RETURN 

J Ni  = N+l-I 

SORT  AVAILABLE  ROhS  BY  DEGREE 

DO  40  J=l, Nl 

J IE  ( J ) = 0 

DO  SO  J- 1 # N 

K 1 = IPR(J) 

1 H ( U ) = K 1 

K2  = NDR(Kl) 

IF  (K2.LE.0)  GO  TO  10 
i IE ( K2  ) = IE ( K2 ) ♦ 1 

DO  bO  J=2rNl 

) IE ( J ) = IE( J)+IE(J-1 ) 

DO  70  J=i*N 

K 1 - I H ( J ) 

K2  = NDR(Kl) 

K3  = 1E(K2)+I-1 
IE ( K2 ) = IE(K2)-1 
) IPR ( K3 ) = K 1 

SORT  AVAILABLE  COLUMNS  BY  DEGREE 

DO  BO  J=1*NI 
I I E ( J ) = 0 

DO  90  J— I i N 
K 1 - IPC(J) 

IH  ( J ) = K 1 
Kid  = NDC(Kl) 

IF  (K2.LE.0)  GO  TO  10 
I IE ( K2 ) = IE ( K2 ) +1 
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00  100  J=2>N1 

100  IE(J)  = IE( J)+I£( J-l) 

00  110  J=1,N 
K 1 = I H ( J ) 

K 2 = NDC ( K 1 ) 

K3  = I E ( K2 ) + I- 1 
IE(K2)  = IE(K2)-1 

110  IPC(K3)  = K 1 
C 

C INITIALISE  FOR  MINIMAL  DEGREE  SEARCH 
C 

1 X - 0 
IOX  = N**2 
JR  = I 
JC  = I 

JRP  = 1PK(JR) 

JCP  = IPL(JC) 

C 

C Tt ST  (-OR  TERMINATION  OF  SEARCH 
C 

1 2 U nDX  = (NOR( JRP)-1)*(NDC(JCP)-1) 

IF  (NDX. Gt .IOX)  GO  TO  240 

IF  (NOC ( JCP) .GT.NDR ( JRP) ) GO  TO  180 


0 73 

C 

074 

C SI  ARCH  IN  THE  COLUMN 

07b 

C 

0 76 

J = JCP 

0 77 

AM  = 0 , 

J 76 

130 

J = CY ( J) 

079 

IF  (J. ED. JCP)  GO  TO  140 

OHO 

AM  : AMAX1 ( AM , ABS ( A ( J ) ) ) 

0 H 1 

oO  TO  130 

0rt2 

140 

AM  = F * AM 

0H3 

160 

J = CY ( J) 

064 

IF  (J. ED. JCP)  GO  TO  170 

06b 

IF  ( ABS ( A ( J ) ) .LT.AM)  GO 

066 

K1  - J 

067 

160 

K1  = R Y ( K 1 ) 

IF  (Kl.Gf.N)  GO  TO  160 
K2  = ( NOR ( K 1 ) - 1 ) * ( NDC ( JCP ) "1 ) 
IF  (K2.GE.I0X)  GO  TO  150 
IX  : J 
KR  =;  K l 
KC  = JCP 
IOX  = K2 
GO  TO  150 
170  JC  = JC  + 1 

IF  (JC.GT.N)  GO  TO  240 
JCP  = IPC{JC) 

GO  10  120 
C 

C StARCH  1 IN  THE  ROW 
C 

160  J = JRP 
AM  = 0. 

190  J = RY ( J) 

IF  (J. ED. JRP)  GO  TO  200 
AM  = AMAX1 ( AM, ABS ( A ( J) ) ) 

GO  TO  190 
200  AM  = F *AM 


J = RYIJ) 

IF  (J, ED. JRP)  GO  TO  230 
IF  <ABS(A< J) ) .LT.AM)  GO  TO  210 


K1  = J 
K 1 = C Y CK1 ) 

IF  (Kl.GT.N)  GO  TO  220 


-^vrr:' 


-■ 
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I 


I 


lib 

117 

118 
119 

iSJ 

122 

123 

124 
12b 
12b 

127 

128 

129 

130 

131 

132 


230 


K2  = (NDR(JRP)-1)*<NDC(K1)-1> 
IF  (K2.bE.IDX)  GO  TO  210 
IX  = J 
KR  = JRP 
KC  = K 1 
1DX  = K2 
bO  TO  210 
JR  = JR  + 1 

IF  (JR.GT.N)  GO  TO  240 
JRP  = IPR(JR) 
bO  TO  120 


C Sb  ARCH  F IN  I bHED » REMOVE  KR  t KC  FROM  AVAILABLE 
C PIVOT  ROWS  AND  COLUMNS 
240  IF  (IX.EQ.O)  bO  TO  10 


DO  2bU  J=I»N 

IF  (IPR(J) .NE.KR)  GO  TO  250 


133 

134 

135 
13b 

137 

138 

139 

140 

mi 

142 

143 

144 

145 
14b 


IF  (J.EO.l)  GO  TO  2b0 
IPR(J)  = I PR ( 1 ) 

IPR ( I ) = KR 

oO  TO  2bU 
250  CONTINUE 
2b0  UO  270  J=1,N 

IF  ( 1 PC ( J ) , NE , KC ) GO  TO  270 
IPC(J)  = IPC(I) 

IPC  ( I ) = KC 
bO  TO  280 
2/0  CONTINUE 
C 

2b0  RETURN 
END 


001 

DU2 

003 

004 

005 

006 

007 

008 
0U9 


C 

C 

C 

C 


SUBROUTINE  NUEC1<M0,FD, A#AN, IE> IH) 

DIMENSION  MD(l),FDU)»A(l)»AN(l)»IE{l),lH(l) 

♦ DECOMPOSE  NONSYMMETRIC  MATRIX  USING  GENERATED  RECORD** 


N = MQ ( 1 ) 
MO (4)  = 0 
FD<5)  = U. 


010 

Oil 

ms 

014 

015 

016 

017 

018 

019 

020 
021 


10 


FD<6) 
FU(7) 
FD  ( 4 ) 
DO  10 
FD  ( 7 ) 
FD(7) 


= 1. 

= 0. 

T 0* 

1=1  fN 

= FD(7)*AL0G<AN<  I ) ) 


= 0.5*Fl(7) 

C CLEAR  EXTRA  STORAGE 

IF  (MD(6) ,LE,MD(b) ) GO  TO  30 
MM  = MQ(b)+l 
Ml  = M D ( 6 ) 

UO  20  I=MM,M1 
2u  A ( I ) = 0, 


I 

L ~ r>#_ 
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024 

02b 

02b 

027 

028 

029 

030 

031 

032 

033 
054 
035 
0 3b 

037 

038 

039 

040 

041 

042 

043 

044 
04b 

046 

047 

048 

049 
ObO 
Obi 

052 

053 

054 

055 

056 

057 

058 

059 
ObO 
Obi 
0b2 

063 

064 
06b 
Obb 
0o7 
068 

069 

070 
0/1 
0 72 

073 

074 
07b 

076 

077 

078 

079 

080 
081 
082 


C INITIALIZE  tl LES 
30  CALL  DRI 
CALL  NVWl 
CALL  NVWd ( - 1 » - 1 ) 

C 

C LOOP  ON  THE  PIVOTS 
C 

DO  130  1 = 1»N 

C GET  PIVOT  ADDRESS  AND  ChECK  PIVOT  MAGNITUDE 
CALL  DR  (IX) 

5 = A( IX) 

IF  (ABS(S)  .LE.FD<1  ) ) GO  TO  150 
FD(b)  = FU(b)+ALOG(ABS(S)  ) 

IF  (S.LT.O. ) FD ( 6 ) = -FD(6) 

FD ( 4 ) = AMAX1 ( FD ( 4 ) # ABS ( S ) ) 

A(  IX)  = 0, 

CALL  DR(KC) 

CALL  DR(KR) 

C Gt T THE  COLUMN  ELEMENTS 
CALL  OR ( K2  ) 

IF  (K2.LE.U)  GO  TO  50 
DO  40  J=1»K2 
CALL  DR(K3) 

6 ALL  URIAH  J)  ) 

A ( J ) = A(*3)/S 

FD ( 4 ) = AM AXl ( FD ( 4 ) # ABS ( A ( J ) ) ) 

4 o A ( K3 ) = 0. 

C OFT  THt  ROW  ELEMENTS 
bO  CALL  DR(K1) 

IF  ( K 1 , Ll  « 0 ) GO  TO  80 
DO  60  J=1*K1 
CALL  DR ( K 3 ) 

CALL  DR( 1H( J) ) 

AN ( J ) = A ( K 3 ) 

FD ( 4 ) = AMAX1 ( FD ( 4 ) #ABS( AN ( J ) ) ) 
bo  A ( K3 ) = 0 • 

C MODIFY  CROSS-POINT  ELEMtiMTS 
IF  (K2.LE.0)  GO  TO  80 
DO  70  J=1»K2 
DO  70  JJ=1»K1 
CALL  DR(K3) 

A(K3)  = A { K3 ) -A  1 J ) * ANtjU ) 

7  0 FD ( 4 ) = AMAX1 (FD(4 ) »ABS( A(K3)  ) > 

C WRITE  OUT  PIVOT  ROW  AND  COLUMN 
8u  CALL  NVWF ( — KC  t -KR ) 

IF  IK2.EO.O)  GO  TO  100 
DO  90  J=1*K2 

90  CALL  NVWF ( IE ( J } » A ( J ) ) 

100  IF  (K1.EU.0)  GO  TO  120 
DO  110  J=1,K1 

110  CALL  NVWri(IH(J)  f AN(JM 

120  CALL  NVWb(-KC»S) 

130  CONTINUE 

c 

CALL  NVWt 
RETURN 
C 

150  MD ( 4 ) = 3 

RETURN 
C 

END 


41  - 


uul 

UU2 

003 

0U4 

OUb 

006 

007 

008 
0 U9 
010 
Oil 
U 12 

01 3 

014 
Ulb 
016 

017 

018 

019 

020 

021 
022 

023 

024 
02b 
026 

027 

028 

029 

030 

031 

032 

033 
0 34 


SUBROUTINE  SSLV ( MD » X # Y ) 
DIMENSION  MD(1)»X(1)»Y(1) 
C 
C 

c 
c 

N = MD{  1 ) 

DO  1U  1 = 1 > N 
lu  r ( I ) = x ( l ) 

CALL  SVRI 
I = 0 

C FORWARD  bACKSUBSTITUTlON 
20  CALL  SVRF ( L2  » Z ) 

I = 1 + 1 
L2  = -L2 
S = Y ( L2 ) 

30  CALL  SVRF(L2»2) 

IF  (L2.LT.0)  GO  TO  40 
Y (L2)  = Y(L2)-Z*S 

00  10  31) 

40  IF  U.LT.N)  00  TO  20 

C BACKwARU  BACKSUBSTITUTlON 
bO  CALL  S vRb ( L2  » Z ) 

1 = i-1 
J = -L2 
Y ( J)  = Y ( J) /l 

oU  CALL  S VHb { L2  t Z ) 

IF  (L2.LT.0)  GO  TO  70 
Y(J)  = Y ( J ) —2  + Y ( L2  ) 
oO  TO  60 

7o  IF  il.Gl.u)  oO  TO  bO 

RETURN 


END 


Oul 

002 

0U3 

004 

00b 

Oub 

007 

008 
009 

011 

012 

013 

014 
01b 
016 

017 

018 

019 

020 
021 


C 

C 

c 

c 


c 

s 

10 

c 

c 

c 


c 

c 


SUBROUTINE  NSLV(MD»X»Yr AN) 

DIMENSION  MD ( 1 ) » X ( 1 ) » Y ( 1 ) » AN ( 1 ) 

* BACKSUBSTITUTlON  FOR  NONSYMMETRIC  DECOMPOSED  MATRIX  * 


EQUIVALENCE  (K5»S) 

N = MD ( 1 ) 

SAVE  RIGHT  SIDE 

UO  10  1=1  * N 
AN ( I ) = XII) 

INITIALIZE  FILES 

CALL  NVRI 
I = 0 
J = 0 

SOLVE  lower  TRIANGULAR  SYSTEM 
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0 22 

C 

023 

20 

CALL  NVKF(K4»S) 

024 

IF  (K4.6T.0)  GO  TO  30 

02b 

I = 1 + 1 

020 

K4  = -K4 

02  7 

K 3 = -Kb 

026 

Y{K4)  = AN ( K3 ) 

029 

IF  ( 1 .Gt.N)  60  TO  40 

030 

01  = AN ( K3 ) 

031 

GO  TO  20 

032 

30 

AN ( K4 ) = AN(K4)-U1*S 

033 

GO  TO  20 

034 

C 

03b 

C 

SOLVE-  UPPER  TRIANGULAR 

0 3b 

C 

037 

40 

CALL  NVRd(K4»S) 

036 

IF  (K4.GT.UJ  GO  TO  SO 

039 

J = J + l 

040 

041 

IF  ( J.Nt. 1 ) Y ( IX)  = Y 
IF  (J.GT.N)  GO  TO  60 

042 

IX  = -K4 

043 

01  = S 

044 

GO  TO  4U 

04b 

bO 

T ( I A ) = Y ( IX )-S*Y (K4) 

046 

C 

GO  TO  40 

047 

048 

bO 

K t T U R N 

049 

C 

ObO 

LNU 

i 

I 


■r- 
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6.2  I/O  Programs: 


* I/O  ROUTINE  FOR  DECOMPOSITION  PROCESSES  * 

******************************************* 

SUBROUTINE  DwI 

THE  DECOMPOSITION  PROCESS  GENERATES  AN  ARRAY  OF 
POSITIVE  INTEGERS.  THIS  PROGRAM  PROVIDES  A BUFFERED 
INPU I /OUTPUT  USING  FILE  10,  THE  ENTRIES  ARE  AS  FOLLOWS: 

DwI  - INITIALIZE  FOR  WRITE 
DW(K)  - WRITE  K AS  NEXT  ENTRY 
DwE  - TERMINATE  WRITING 
OR!  - INITIALIZE  FOR  READ 
OR ( K ) - READ  NEXT  ENTRY  K 

IX  IS  THE  BUFFER  SIZE 
PARAMETER  lx  = 2S0 
DIMENSION  IB (IX) 


+ initialize  for  writing  * 


j = l 
REWIND  1U 
RETURN 

ENTRY  OW(K) 

* WRITE  NEXT  ENTRY  K * 


I B ( J ) = K 
J = J+l 

IF  (J.LE.IX)  RETURN 
WRITE  (10)  It) 

J = 1 
RE  TURN 

ENTRY  OWE 

* TERMINATE  WRITING  * 

********************* 

IF  (J.NE.l)  WRITE  (10)  IB 
RETURN 

ENTRY  DRI 

********************** 

* INITIALIZE  READ-IN  * 
********************** 

REWIND  10 
J = IX 
RETURN 

ENTRY  DR ( K ) 

* READ  NEXT  LNTRY  K * 

********************* 

J = J+l 

IF  (J.LE.IX)  GO  TO  10 
READ  (10)  IU 
J = 1 

i K = IB ( J ) 

RETURN 


T" 


w 
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uu  l 

002 

003 

004 
00b 
006 
007 
OOd 

009 

010 

011 

012 

013 

014 
01b 
Olo 
017 
Old 

019 

020 
021 
022 

023 

024 
02b 
026 

027 

028 

029 

030 

031 

032 

033 


036 

037 
03d 

039 

040 

041 

042 

043 

044 
04b 

046 

047 

048 

049 
ObO 
Obi 
0b2 
0b3 
Ob4 
Obb 
Ob6 
0b7 
0b8 
Ob9 
060 
061 
062 
063 
Ob4 
06b 
Ob6 

067 

068 

069 

070 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

L 

C 

c 

i 

c 

c 

c 

c 

0 


c 

c 

c 

c 

o 

0 

0 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


*********************************************** 

* I/O  HOUTINE  FOR  DECOMPOSED  SYMMETRIC  MATRIX  * 

*********************************************** 


SUBROUTINE  sywl 
DECOMPOSED  MATRIX 


lb  PLACED  IN  FILE  11  AS  A RANDOM 
TS  OF  A DOUBLE  ARRAY  WHICH  IS 


THE 

ACCESS  file,  it  consi 
BUFFERED.  ALTHOUGH  THE  FIRST  PART  OF  THE  ARRAY 
IS  AN  iNitGER  (SIGNED)  ARRAY » THIS  ROUTINE  DOES  NOT 
PACK  IT.  THE  ROUTINE  HAS  THE  FOLLOWING  ENTRIES: 


S vw  I - INITIALIZE  FOR  WRITE 

S V W ( A 1 » A2  ) - WRITE  A1.A2  AS  NEXT  ENTRY 

SVwE  - TERMINATE  WRITE 

SVR1  - INITIALIZE  FOR  READ 

SVRF ( A1 » A2 ) - READ  NEXT  ENTRY  A1#A2 

SVRd ( A1 » A2 ) - READ  PREVIOUS  ENTRY  Al»A2 


ThE  ROUTINE  ASSUMES  THAT  THE  WRITTEN  ARRAY  Is  READ  ONCE 
FORWARD  THEN  READ  BACKWARD. 


PARAMETER  NX  = 100 
PARAMETER  MX  = 280 
PARAMETER  MXX  = 2*MX 
NX  IS  ThE  MAXIMUM  NUMBER  OF  RECORDS 
MXX  IS  THE  LENGTH  OF  THE  RECORDS 
DIMENSION  B ( 2 » MX  ) 


* INITIALIZE  WRITING  * 


N = 0 

J = 1 

DEFINE  FILE  11 (NX#MXX,U#  IX) 

IX  = lx 

return 

ENTRY  SVW(A1»A2) 

* WRITE  NEXT  ENTRY  A1.A2  * 
************************** 

8 ( 1 » J ) = A 1 
d ( 2 » J ) = A2 
J = J+l 

IF  (J.LE.MX)  return 
N = N+l 

WRITE  (11 'N)  V 

J = 1 

RETURN 

ENTRY  SVWt 

********************* 

* TERMINATE  WRITING  * 
********************* 

IF  (J.EU.l)  RETURN 

N = N+l 

WRITE  < 1 1 * N ) B 
RETURN 

ENTRY  SVRI 

********************** 

* INITIALIZE  REAO-IN  * 
********************** 

M = 0 
J = MX 
RETURN 
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071 

C 

072 

ENTRY  SVHf-  ( A1  , A2) 

073 

c ************************* 

074 

C * REAL)  NEXT  tNTRY  A1*A2  * 

07b 

c ************************* 

0 76 

C 

077 

J = J + l 

0 78 

IF  (J.LE.MX) 

GO  TO  10 

0 79 

M = M+l 

060 

KEAU  (ll'to) 

B 

Obi 

J = 1 

0H2 

10  A 1 = H(lrJ) 

063 

A2  = B ( 2 » J ) 

084 

ke  turn 

065 

c 

066 

entry  SVRB(A1,A2) 

067 

c ***************************** 

068 

U * READ  PREVIOUS 

ENTRY  A1,A2  * 

089 

C ***************************** 

090 

c 

091 

IF  (J.GT.O) 

GO  TO  20 

092 

to  = M— i 

093 

READ  (11 'to) 

B 

094 

J = MX 

095 

2U  A 1 = B ( 1 mJ  } 

096 

A2  = 6 ( 2 # J ) 

097 

J = J-l 

ova 

RETURN 

099 

C 

100 

END 
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UU1 

002 

0u3 
004 
ddb 
00b 
007 
Odd 
009 
did 
dll 
d 12 
d 1 3 
d 14 
d lb 
d lb 
d 1 7 
did 

019 

020 
021 
022 

023 

024 
02b 
02b 
027 
0 2d 

029 

030 

031 

032 

033 

034 
03b 
0 3b 
037 
03d 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

u 

c 

c 

c 

c 

c 


♦ I/O  ROUTINE  FOR  NONSYMMETRIC  DECOMPOSED  MATRIX  * 

>****t ********** 

subroutine  nv w i 

I HE  LOWER  AND  UPPER  TRIANGULAR  MATRICES  OF  ThE 
Ufc  COMPOSED  NONSYMMETRIC  MATRIX  ARE  CONTAINED  IN 
FILE  12  ANO  FILE  13*  RESPECTIVELY,  AS  RANDOM 
ACCESS  FILES.  THEy  ARE  IN  THE  FORM  OF  BUFFERED 
DOUBLE  ARRAYS.  THE  ENTRIES  ARE  AS  FOLLOtoS, 


NVwI  - INITIALIZE  FOR  WRITE 

N V toF  ( A 1 # A2 ) - WRITE  A1,A2  AS  NEXT  ENTRY  ON  FILE  12 

NVWb(Al,A2)  - WRITE  A1,A2  AS  NEXT  ENTRY  ON  FILE  13 

NVwE  - TERMINATE  WRITING 

NVRi  - INITIALIZE  FOR  READ 

N vRF ( A i » A2 ) - READ  NEXT  ENTRY  FROM  FILE  12 

NVRB ( A 1 , A2 ) - READ  PREVIOUS  ENTRY  FROM  FIlF  13 

FILE  12  IS  READ  FORWARD,  FILE  13  BACKWARD, 

PARAMETER  NX  = 100 
PARAMETER  MX  = 2d0 
PARAMETER  MXX  = 2*Mx 
NX  IS  THE  MAXIMUM  NUMBER  OF  RECORDS, 

MXX  IS  THE  RECORD  SIZE 

DIMENSION  B12(2,MX ) ,B13(2,MX) 


* INITIALIZE  WRITING  * 


N12  = 0 
N13  = U 
J12  = 1 
Ji 3 - 1 

DEFINE  FILE  12(NX, MXX, U, 1X12) 
1X12  = 1X12 

DEFINE  FILE  13(Nx, MXX, U, 1x13) 


0 09 

040 

041 

042 

043 

044 
04b 
04b 

047 

048 

049 
ObO 
Obi 
0b2 
0S3 
Ob4 
Obb 
0b6 
0b7 
Obd 
Ob9 
dbd 
061 
062 
0b3 

064 

Obb 

066 

0b7 

Obd 

0h9 


1X13  = 1X13 
RETURN 
C 

ENTRY  NVWF(A1,A2) 

C * WRITE  A 1 , A2  ON  FILE  12  * 

c ************************** 

c 

B12 ( 1 , J12 ) = A 1 
B12(2,J12)  = A2 
D12  = J12+1 
IF  (J12.LE.MX)  RETURN 
in|  12  = N12  + I 
J12  = 1 

WRITE  ( 12 ' N12 ) B12 
RETURN 
C 

ENTRY  NVWB ( A1 , A2 ) 

C ************************** 
C * WRITE  A1,A2  ON  FILE  13  * 
C ************************** 

C 

B13(1,J13)  = A 1 
dl 3 ( 2 , J13 ) = A2 
J13  = J13+1 
IF  (J13.LE.MX)  RETURN 
N13  = N13+1 
J13  = 1 

WRITE  (13*N13)  Bl3 
RETURN 
C 
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U 70 
U /I 
072 
U 73 
U 74 
0 7b 
U 76 
J 77 

078 

079 

oao 

Odl 

082 

063 

084 

08b 


08t> 

Ob7 
Odd 
od9 
0 90 

091 

092 

093 

094 
09b 

096 

097 
09d 


099 


100 

1U1 

102 

1 03 

104 


lub 
A 06 
1 0 7 
108 

109 

110 
111 
112 

113 

114 


ENTRY  NVwE 

C ********************* 

C * TEKMINATE  WRITING  * 

C ********************* 

C 

IF  (J12.E0.1)  GO  TO  10 
N12  = N12  + 1 
WRITE  (12'N12)  B12 
10  JJ  = Mx 

IF  (J13.EQ.1)  RETURN 
N13  = N13  + 1 
WRITE  ( 1 3 1 N1 3 ) B13 
J J — J 1 3” 1 
RETURN 
L 

cNTK  Y NVRi 

C ********************** 

C * INITIALIZE  READ-IN  * 
c ********************** 


N2  = 0 

N3  = N13+1 
J2  = MX 
J3  = 0 
RETURN 

ENTRY  NVRF ( A1 » A2 ) 

.***************** 

* READ  A 1 * A2  f-ROM  FIlE  12  * 


2 0 


C 

c 

c 

c 


J2  = 02+1 
IF  (J2.Lc.MX) 
N2  = N2  + 1 
READ  <12'N2) 

J2  = 1 

A1  = B 12 ( 1 » J2  ) 
A2  = Bl2(2fJ2) 


GO  TO  20 
B12 


RETURN 

ENTRY  NVRB(A1,A2) 

******* 

♦ READ  PREVIOUS  ENTRY  Al,A2  FROM  FILE  13  * 
****************************************** 


lib 

116 

117 

lid 

119 

120 
121 
122 

123 

124 

125 


3u 

C 


IF  (J3.GT.0)  GO  TO  30 

N3  = N3-1 

READ  ( 1 3 * N 3 ) Bl3 

J3  = MX 

IF  (iN3.E0.N13)  J3  = JJ 
A 1 = B13 ( 1 » J3 ) 

A2  = B13 ( 2 » J3  ) 

J3  - J3-1 
RETURN 

end 
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